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Abstract 

We investigate the density matrix renormalization group (DMRG) discovered 
by White and show that in the case where the renormalization eventually con- 
verges to a fixed point the DMRG ground state can be simply written as a 

"matrix product" form. This ground state can also be rederived through a 
simple variational ansatz making no reference to the DMRG construction. We 
also show how to construct the "matrix product" states and how to calcu- 
late their properties, including the excitation spectrum. This paper provides 

details of many results announced in an earlier letter. 
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I. INTRODUCTION 



After Wilson's development of the renormalization group (RG) to solve the Kondo 
problemS it was believed that RG could be used for other problems as well. Kadanoff's 
blocking technique combined with Wilson's RG idea was applied to problems like quantum 
lattice systems such as the Hubbard and Heisenberg models but progress turned out to be 
surprisingly difficult. However, in 1992 White developed the density matrix renormalization 
group (DMRG)iH method which since then has had a spectacular success in calculating 
ground state energies and other static properties of many ID quantum systems. In this 
paper we explore the nature and underlying principles of the DMRG to find out why the 
results of DMRG calculations are so remarkable accurate. A summary of this work has been 
presented in an earlier paper!, and the present paper provides a complete discussion and 
derivations of the results. For background information on the DMRG there are excellent 
articles by White!!. 

In Sec. H we give a very brief summary of the DMRG. In Sec. |Hj we show that if the 



DMRG algorithm converges to a fixed point, the DMRG ground state leads to a special 
ansatz form for the wave function, demonstrating the equivalence of the DMRG to a varia- 
tional calculation. To make things more concrete we apply our ideas to the antiferromagnetic 
Heisenberg spin-1 chain with quadratic and biquadratic interactions, defined by 

J fJ = ^S J -S m -/3(S l -S i+1 ) 2 . (1) 

i=i 

In Sec. [IV] we define a set of variational states using the special ansatz form of Sec. |T|. 
In Sec. [V| we extend the ansatz to include a set of Bloch states that describe elementary 
excitations in both finite and infinite systems. These calculations are rather lengthy and the 
details can be found in appendices. Sec. [VJ contains some numerical results for the spin-1 
chain to compare our variational ansatz to more involved calculations. 

We would like to mention that all numerical work described here was programmed with 
Mathematica on an ordinary desktop workstation. Each calculation we describe here take 
anywhere from a few seconds to a few minutes. 



II. DENSITY MATRIX RENORMALIZATION GROUP (DMRG) 

Since the DMRG was discovered by S.R. White! in 1992 it has had great success in 
describing ID interacting quantum systems!'! - !. Ground state and excited state properties 
have been calculated to high accuracy with modest computational effort. With hindsight, 
it will be seen that the ideas of this paper do not logically depend on the DMRG, but they 
were inspired by the DMRG and we will therefore begin this section by summarizing some 
aspects of the DMRG. 

In a renormalization scheme like the DMRG one typically starts with a very short ID 
chain and then lets the length increase by iteratively adding a single site. After each new 
site, an approximate Hamiltonian is constructed. This is done by keeping only a small 
subspace of the Hilbert space to keep the Hilbert space at a manageable size as one lets 
the chain grow. The central idea in DMRG is to keep the "most probable" states when 
truncating the basis in contrast to the usual old-fashioned real space RG methods (see e.g. 
Ref. ^ and references therein) where the lowest energy eigenstates are kept. The way to 
achieve this is to split a complete system ("universe") into two parts, a "subsystem" and 
an "environment" and then to construct the reduced density matrix for the "subsystem" as 
part of the "universe" . The state of the "subsystem" is then given by a linear combination 
of the eigenstates of the density matrix with weights given by the eigenvalues. 

The renormalization starts with a short ID lattice with just a few sites. Label this system 
Hb- A renormalization step of the DMRG can be described by the following algorithm: 

1. Construct the Hamiltonian for the "universe", H$ = Hb + Hi + H^ + H§, where 
Hb comes from the previous iteration and Hi is a new site added. The superscript R 
denotes a second block that is reflected before joined to the other parts. The block Hb+ 
H\ now is our "subsystem" and H^ + H§ our "environment". The Hamiltonian matrix 
Hs for the "universe" is constructed with tensor products involving the intrablock parts 
Hb and Hi and the interactions between the blocks. 
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2. Diagonalize H s to obtain the ground state of the "universe". This state is called 
the target state. 

3. Construct the reduced density matrix p iiV = Y.j where |^) = ® \j) 
and \ j) are basis states of the "subsystem" and the "environment" respectively. 
The eigenstates of p with the highest eigenvalues correspond to the most probable 
states of the "subsystem" when the "universe" is in the state |\&). 

4. Now choose the m states of the diagonalized density matrix with highest eigenvalues to 
form a new reduced basis for the block H B + Hi . Project the Hamiltonian and other 
operators onto this basis by Hb> = A{Hb+i)A\ where A is the projection operator one 
constructs from the kept eigenstates of the density matrix from step 3 and H B +\ is the 
Hamiltonian matrix for the "subsystem" . If the single site added has m s states in its 
basis, A is represented by a m x (m * m s ) matrix and Hb+i by a (m * m s ) x (m * m s ) 
matrix. 

5. Rename H B ' to H B . 
This completes an iteration. 

III. THE MATRIX PRODUCT STATE 

To begin the renormalization procedure one starts with a block consisting of a short lat- 
tice whose basis states can be calculated exactly. When the renormalization proceeds and 
the chain described by the block gets longer we don't use the full set of basis states for de- 
scribing the block but have to discard some part of the Hilbert space in each renormalization 
step. 

Assume we have a block that represents a chain with n — 1 sites. Let m s be the number 
of possible states of a singe lattice site. If we would treat this system exactly there would 
be m™" 1 states in the Hilbert space basis for this system. In the case of a spin 1 chain, we 
could label the site with the z-component of the single spin-1, so that m s = 3. The number 
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of states in this complete basis rapidly becomes too large to handle when n is increased. 
Assume therefore that an approximation is made and our chain is represented by a smaller 
set of states labeled by This set of states has been chosen by the previous iterations 

of the renormalization with the aim to describe the low energy physics. Assume there are m 
states in this basis, where m < m" -1 . If this is the first iteration, {\j3) n _i\ is the complete 
basis. 

We now add a single site, labeled by s n , the ^-component of spin, to the left hand side 
of our block resulting in a new block with n sites and m s x m states in its basis. The 
basis states are now generated by the product representation {\s n ) <g> |/3) n _x}. We now use 
a projection operator A n to generate a new truncated basis with typically m states that 
represent the "important" states of the longer block. This whole process is written (see also 
Fig-D 

\a) n =J2K' iM K)®\P)n-i, (2) 

where we have indexed A by the chain length n and its matrix indices a and (/?, s n ). Note 
that (/?, s n ) is thought of as a single index labeling a tensor product of the states \s n ) and 

| 

In the DMRG, a specific algorithm is used to calculate A, but this is not important in 
the present discussion. 

We now make two crucial observations. 

1. First we perform a simple change in notation: A" ,|8 [s n ] = A"'^' Sn \ thus writing the 
m x (m s * m) matrix as a set of m s m x m matrices. 

2. Second, we assume that the recursion leads to a fixed point for the projection operator 
so that we can write A n [s] — > A[s], as n — *■ oo. 

By recursively applying the renormalization step in Eq. (0) we now find that 

\a) n = £ (A[s n ]A[s n . 1 ]...A[s 1 ]) a ^\s n s n . 1 ...s 1 }®\(3) , (3) 

s„...si,{3 
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where |/3)q represents some state far away from n. We thus see that the renormalization 
procedure results in a wave function that can be written in a matrix product form. Eq. (|3|) 
now suggests a natural form for the wavefunction with the following ansatz. 
For every m x m matrix Q, we define the (unnormalized) state \Q) n 

\Q) n = £>( QA[s n ] . . . A[s x ] )\s n ...s x ). (4) 

M 

Thus \Q) n can be viewed as a state that is uniform in the bulk, but with a linear combination 
of boundary conditions defined by \a) n on the left and \/3) on the rightH. The special case 
of Q = 1, the identity matrix, leads to a state with periodic boundary conditions. This 
Q = A. state we will later on use as our trial ground state. 

If we now demand that the projection of Eq. (0) preserves orthonormal bases, {a\a') = 
5 a , a ', we can use the recursion formula Eq. (fj) and the orthogonality of the local spin states 
and previous block states to find 

<W= E (A a '^[s']YA a ^[s} (s'\s)(P'\/3) 

/3,0',s,s' 

= Y,(A{s]As]r a '- (5) 

s 

Hence in matrix form we have J2 S A[s]A*[s] = 11. This constraint will be used later to reduce 
the number of free parameters in A. 

We now analyze the projection matrix A. The Hamiltonian of Eq. (HD is spin rotationally 
invariant since it commutes with all three components of the total spin S to t — J2i Si- In order 
that the projection in each step preserves this symmetry, our basis states of a block must 
form a representation of total spin. Since we keep basis states with many different values of 
total spin as well as many states with the same total spin in each iteration, all the basis states 
together must form a sum of irreducible representations of total spin. Adding a spin one 
does not mix even or half-odd spin representations, thus the basis states must form a sum of 
either all half-odd or all integer spin representations. Most naturally for the spin-1 chain one 
would work with integer spin representations, but by placing a single spin-1/2 on the right 
hand side of the entire chain one could use half integer spin representations to represent the 
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blocks instead. This is consistent with the existence of a spin-1/2 edge stateua and we have 
found that working with half-odd integer representations give far better numerical results. 

We now discuss how the projection operator A can be constructed. In our numerical 
work, we have kept 12 basis states in each iteration and we have used the half-odd spin 
representations. By doing a DMRG calculation on the spin-1 chain we have found that when 
approximately 12 states are kept the blocks are represented by a sum of two spin-1/2 and two 
spin-3/2 irreducible representations. Since there are two sets of each representation, we have 
to introduce a new label, 7, to distinguish them. The "old" representations representing the 
old block we have uniquely labeled by ordinal number 7, with the corresponding total spin 
j. (See Fig. ^.) Implicit in the labeling of the states is the z-component of total spin m, 
which is not shown in Fig. 0. These are thus the twelve "old" states I7, m), that fall into 
the four different irreducible representations of total spin. 

After adding a single site and then truncating the Hilbert space we get the "new" basis 
states similarly labeled by 7' and their corresponding total spin j', to the right in Fig. |2|. 
These will thus be the twelve basis states {'y' ,m') that represent the new block with one 
more site. 

Let us now examine what happens in our example when going from the old 7 to the new 
7'. When adding a single spin-1 to the old block of twelve states we get 36 "intermediate" 
states in the product representation of the old block states with a spin 1. These states fall 
into 10 irreducible representations labeled in Fig. Q by their total spin j". We then project 
from these 10 reps back down to the four reps that we have chosen to keep. This projection 
must preserve the total spin symmetry, i.e. it can't mix different j" and cannot depend on 
total m". The only nonzero projection terms P 7 '> 7 are indicated by lines in the figure. Since 
there is exactly one "intermediate" spin-1/2 and one spin-3/2 for each of the four "old" 
representations 7 there is one projection term from each of the "old" 7 to each of the "new" 
7'. There are thus 16 nonzero projection terms which are in fact not independent, but are 
related by the requirement that the new states are orthonormal. 

Expressing all this mathematically, we let, as in the figure, 7 uniquely label a rep of total 
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spin of a block and 7(7) denote the value of total spin of that rep. Each state is thus labeled 
by 1 7, m) where m is the z-component of total spin. The single spin to be added is labeled 
by \s), where s is the ^-component of the spin-1. The new states are thus given by 

IYX> = E^"l7,j(VW>, (6) 

where I7, .7(7')) m') denotes the 36 intermediate states formed by \s) ® I7, m) written in 
the total spin basis. Since we demand that the projection preserves total j and m, these 
states can be explicitly constructed using the Clebsch-Gordan coefficients on the form 

((Ju m i)U2,m2)\j,m) as 

W,j(j'),m') = Y,(U(l),rn)(l,s)\j(Y),m>) (\s) ® | 7 » . 

m,s 

Inserting this into Eq. (§) we find that 

\i,m')= £ A (l '' m '^^[s](\s)® | 7 ,m», 

s,(7,m) 

where 

A (y,m'),(7,™) W =P 7',7 ((j(7);m) (i, s )| j( 7'),m'). 

Thus, although the projection matrices A contain a total of 3 x 12 x 12 numbers, they are 
in fact generated by the relatively few degrees of freedom available in P 7 '> 7 . 

For this case with twelve basis states there are naively 16 parameters in P 7 ,7 . Demanding 
normalization of all basis states, (71, mil 72, m 2 ) = 5 7l , 72 5 mi , m2 , yields the condition that the 
diagonal elements of P T P are all 1, where the superscript T denotes transpose. This gives 
four constraints. (Cf. Eq. (^).) Then the basis states of the two spin- 1/2 and the two 
spin-3/2 must be orthogonal, yielding the condition (P T P)Ti'f 2 = 0, whenever j'(7i) = J (72) 
with 71 7^ 7 2 . This gives two more constraints. The spin-1/2 basis states are automatically 
orthogonal to the spin-3/2 states. Finally, a unitary transformation can mix the two spin-1/2 
and likewise the two spin-3/2. Without loss of generality we can fix this freedom, yielding 
two more constraints. We thus end up with only eight free parameters^ in P 7 '' 7 . In the 
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simpler case of saving only six basis states, only two free parameters are available by similar 
arguments. 

With only a few free parameters we can use a variational principle for the energy to 
determine these. At this point it is clear that the DMRG plays no essential role in the 
construction aside from providing a guide to which representations to keep. Even this 
choice could be done variationally. 

IV. THE SET OF STATES \Q) 

A. The ground state ansatz 

To do the variational calculation we need an expression for the energy. As an ansatz 
for the ground state wave function we take the translationally invariant state Q = 11 from 
Eq. (|J) which we denote by Thus 

\l) = J2^(A[s n ]...A[s 1 ])\s n ...s 1 ). (7) 

{*j} 

Note that although it is not explicitly written out, |1) has a definite number of lattice sites 
n. We note that (1|1) = 1 due to Eq. (|5]). For the AKLT modeli ((3 = —1/3) our ground 
state ansatz is exact as are the "matrix product" states of Ref. p7T| [T~5|. 

The expectation value of an operator h, e.g. energy or correlation function, in this state 
is given by 

(l\h\l) = J2 tr ( A *M ■ ■ ■ A*[ Sl '] ) tr ( A[s n ] . . . A[ Sl ] )(s n ' . . . Sl '\h\s n ... Sl ). (8) 

To write this expression in a simpler form we define the tensor product matrix (B<S>C) by 
(B (g) Cp a ^'^ T ' v ^ = B a ' T C l3 ' u . We will in the rest of the article interchangeably use ordinary 
matrix indices a, (3 and composite indices (a, (3), where composite indices are written with 
a parenthesis around them. This means that we can write a m x m "matrix" A as either 
a matrix A a ^ or as a m 2 vector A^ a ^\ When the indices are not explicitly written out, 
the matrix or vector character of the symbol is assumed to be clear from the context. 
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We now use the trace and matrix product identities tr ( B ) tr ( C ) = tr ( B ® C ) and 
(BCD) <g> (EFG) = (B®E)(C® F)(D ® G) to find 

= £ tr ( ( A * W ® A[s n ]) . . . ® A[s x ]) )(s n ' . . . si>|s n . . . si). (9) 

{'iU'i'} 

To write this in a more compact form we define a mapping M from 3x3 spin matrices 
M to m 2 x m 2 matrices M by 

M = J2M s , jS (A*[s']®A[s}). (10) 

We denote by S = (S^, 5 y , S 12 ) the spin-1 representation of total spin and thus by S = 
(S x : S y ,S z ) the "hat" mapping of the 3x3 spin matrices 5". By 1 we denote the "hat" 
mapping of the 3x3 identity matrix. We now see from Eq. (|9]) that the norm and the 
expectation value of the spin at the site j is given by 

(1|1>= tr(l") 
(1|S,|1>= tr(l«-\S), 

where we in the last equation have used the cyclicity of the trace. Other expectation values 
are also easily obtained. Since we can factorize matrix elements like 

s jf s j|Sj ■ Sj|sj,Sj) = (Sj • Sj) s '^ jS / iiS ^. jSi 

we find that expectation values of energy and spin-spin correlation function are given by 

(1|Sj--Sj + i|1) = tr( l n - 2 SS) 

(i\Sj ■ s j+l \i) = tr ( r-^sV^s). (11) 

Similar formulas have also been derived by Fannes ei.a@. 

A more complicated operator, like the biquadratic term (Sj • S.,) 2 , does not factorize 
as neatly and we cannot write the expectation value in such a nice form as above. For 
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these cases we have to replace the term SI 1 1 S inside the trace in Eq. (|Tl|) by the more 
complicated expression 

£ (af j} s'^Si • S,) 2 |s„ Si ) (A[s'i] ® AisiVV-^Atfj] <8> Ay ). (12) 

In order not to make the equations unreadable by crowding them with indices we will in the 
rest of this chapter only present formulas for the ordinary bilinear Heisenberg Hamiltonian 
{(3 = in Eq. ([]])). An interested reader can then generalize the formulas to include the 
biquadratic term, without any fundamental difficulties. 

An important quantity is the string correlation functionB defined by 

9® = (Ss(u^)sf). (13) 

Although the spin-1 chain does not have long range Neel order, it is believed to have a hidden 
long range order that is characterized by the string correlation function. In our ground state 
|1) it is easy to show that it is given by 

g{l) = tr ( P'^S^e^) 1 ' 1 ^ ). 

We note that the spectrum of correlation lengths, i.e. the collection of all possible ex- 
ponential decay lengths £ of correlation functions of the form (Oi{x)0>z{y)) oc e~' x ~ y ^^, is 
determined by the eigenvalues of 1. One can show that 1 is guaranteed to have an eigenvalue 
of 1 due to Eq. (^j), and numerically we find that all other eigenvalues have absolute value 
less than 1. It is however not true that the eigenvalue 1 will always dominate. If each of the 
rows of 0\ or each of the columns of O2 is orthogonal to this particular eigenvector, another 
eigenvalue will determine the correlation length. Thus, the correlation length £ is given by 

t = -r-> ( 14 ) 

logx 

where x is the largest eigenvalue of 1 not orthogonal to the operator. Since the rows and 
columns of the spin operator S turn out to be orthogonal to the eigenvalue 1 while the 
next leading eigenvalue is not, the next largest eigenvalue will determine the decay of spin 
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correlations. The string operator e t7TS i of Eq. fll3| ) turns out to have the same eigenvalue 
spectrum as 1. This time however, the eigenvalue 1 of J is not orthogonal to e l7TSz , giving 
the long range string correlations. 

A possible problem with the construction of the projection operator is that parity is not 
built into the construction of the ground state since the projectors operate from the left to 
the right. There is therefore the possibility that parity is violated in the ground state |1). 
We now investigate this possibility and show how parity is maintained. 

Let V be the parity operator. We thus have 

V\\) = tr(A[s n ]...A[ Sl ] )V\s n ... Sl ) 

M 

= £ tT(A[s n }...A[ Sl } )\ Sl ...s n ). 

M 

Suppose now that there exists an invertible m x m matrix Q-p such that 

Q v A[s] = Si gn[V](A[s]) T Q P , (15) 

where A T denotes transpose and sign[V] is a proportionality constant that will be seen to 
be the eigenvalue of the parity operator. Then it follows that 



P\l) = tr ( QvQvA[s n ] . . . A[ 8l ] )\ Sl ...s n ) 

{*j} 

= sign[V} n tr ( AT [sn] ■ ■ ■ A T [ Sl ] ) |s a . . . s n ) 

= sign[V\ n Y, tr(A[ Sl }...A[s n ] )\ Sl ...s n ) 

M 

= sign[P] n \l). 

Thus, for the ground state to have definite parity, it is sufficient that such a Q-p exists. How 
do we find this matrix, if it exists? We multiply both sides of the defining relation Eq. fllBl) 
by A^[s] and sum over s. Using Eq. (|) we find that 

s 

= sign[P] fe^r^'ls]) QV 
= sign[V) E(^ T W ® A[s])^^Q T /. 

s 
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Thus, Q-p, if it exists, is the eigenvector of the matrix Zls(^4 T [ s ] <8>A[s]) with eigenvalue ±1. 
It is straightforward to numerically obtain the eigenvalue spectrum of this operator, and in 
the cases that we have looked at, this parity operator exists. 



B. The general state \Q) 

We now analyze the set of states \Q) n , for general Q. These states can be interpreted as 
states homogeneous in the bulk but with nonuniformity near the boundary. 

To calculate the norm we use the same trace and tensor product identities as when 
deriving Eq. (|9[). We find that 

(Q'\Q)n = E tr ( (Q'YA*[s n ] . . . A*[ Sl ] ) tr ( QA[s n ] . . . A[ Sl ) ) 
{*,} 

= tr( (Q'*®Q) 1" ). (16) 



We can rewrite this trace as ordinary matrix products. To do this we first define the 
generalized transpose M Tp i- p 2.P3>P4 of a matrix M by 

where {pi,P2,P3,Pa} is a permutation of {1,2,3,4}. We also define a tilde operator M by 
the formula 

tr((Q / ®g)M)= E (Q , ) a '' /3 '(M) (a '' /3 ' ) ' (a * /3) Q a ' /3 , (18) 

a',P',a,P 

so that the tilde operator effectively generates the matrix corresponding to the inner product 
of Q' and Q with M. One finds by writing out Eq. flTJD in components that 

M = M Tsi42 . 

Hence 

((QT''"'\Q a '% = Q' {a '' p,) G(nY a '^ a ^Q^\ (19) 

with 
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G(n) = (H (20) 

The nice thing about Eq. ( |i~9"|) and (p0|) is that we have effectively turned the computation 
of the trace in Eq. ([16]) for all Q and Q' into a matrix inner product between Q, Q' and a 
single m 2 x m 2 matrix G, independent of Q and Q' . Note that on the right side in Eq. ( |i~9"D 
we write Q and Q' as vectors of length m 2 . 

Similarly we can compute the expectation value of the Heisenberg Hamiltonian defined 
in Eq. (0) with (3 = as 

ra-2 

(Q'|#op|Q)n = E tr ( (Q' ® ), (21) 

where 5* denotes the hat mapping in Eq. flTU| ) of the spin-1 matrices. The z-component of 
total spin, (S^) op = J2i Sf is given by 

n-l 

(Q'\(S^) op \Q) n = tr ( (Q' ® Q^l"- 1 " 1 ). (22) 

If we have a more complicated Hamiltonian, like Eq. ([!]) with (3^0, the Hamiltonian 
matrices SS in Eq. (|21|) must be replaced by an expression similar to Eq. fll2|). As we did 
with the norm in Eq. (|l^) we can rewrite Eq. (pip and ( p^) as matrix products by putting 
the summations inside the traces and by using the tilde transformation of Eq. (|IBD, yielding 

(Q'\H op \Q) n = Q'H{n)Q (23) 
(Q'l(^)o P |g)n = Q'^Hg, (24) 

where 

ra-2 

jf (n) = E (i^Sl"- 2 ^) (25) 

i=0 
ra-1 

S£(n) = £ (i^i"" 2 ^) , (26) 

8=0 

where the tilde symbols indicate that the transformation in Eq. fll8|) should be performed 
on the whole sum. 
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In Eq. ( |19~D we determined the expression (Q'\Q) n = Q'G(n)Q for the inner product of 
the states \Q) n and in Eq. (f23|) we found (Q'\H op \Q) = Q'H{n)Q. Since G turns out not 
to be proportional to the identity matrix, we see that the naive basis states, i.e. the states 
(Qi,j) a '^ = 8i,a$j,i3, with i = 1 . . . m and j = 1 . . . m, are not orthonormal. It is not only 
convenient to have an orthonormal set of states, we also want them to be eigenstates of the 
Hamiltonian. The energy of the state \Q) n defined by the m x m matrix Q is given by 

_ QH{n)Q 

" QGin}Q- (27) 
The eigenvalue equation we have to solve is thus 

H(n)Q = E Q G(n)Q. (28) 

We will now construct a set of states that are orthonormal and satisfies Eq. Since 
G is Hermitian we can define a unitary matrix V by the transformation that diagonalizes 
G: 

V^GV = D G , (29) 
where Dq is a diagonal matrix. We now define 

u = V(D G )- 1 ' 2 

so that v)Gu = 1, the identity matrix. We also define 

h = v)Hu (30) 
s z = u^S^u, (31) 

with H and from Eq. fl25|) and Eq. (|26|). It can be verified that [h, s z ] = so that both 
total spin and the energy can be diagonalized simultaneously. Numerically we diagonalize 
h + es z where e is a small number, so that w^(h + es z )w = E + es z is diagonal and we find 
that both h and s z are thereby diagonalized by 

w^hw = E (32) 
w^SzW = m z , (33) 
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with w'w = 11 and where E and m z are diagonal matrices containing the energy eigenvalues 
and the eigenvalues of total spin respectively. Putting Eq. ( PDD and into Eq. (|32|) and 
(|33|) we see that 

(( u w) j ) , H a ^(uw)p n = £ 7 5y, 7 (34) 

a,/3 7 ,a 

Z) (( MW ) t ) , Q G a,p{ UW )p,1 = <V, 7 , (35) 

where a, /3, 7 and 7' are matrix indices. Thus, the columns of (uw) contain the orthonormal 
eigenvectors of H{n) and S^{n). Combining Eq. (j34"D and (|35|) we find 

t p{uw)p n — E y G a) p(uw)p >T 

Hence the matrices 

{Q-y)a,p = W(«,(3) )7 , (36) 

where Q 1 are m 2 m x m matrices, are orthogonal with respect to G and are simultaneous 
eigenstates of H and S^. We therefore define the orthonormal set of states I7) we were 
looking for by 

|7> = IQt)- 

To summarize, we finally have 

(7' I H op I 7) = Ey8yl n 

(TK^^lT) = (^) 7 <V,7 

(t'It) = <W 

The states I7) form a natural basis for describing edge states in finite size calculations, a 
feature which is not further explored in this manuscript. Nevertheless, we will benefit from 
this derivation in the next section were a set of Bloch states are defined in a similar manner. 
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V. BLOCH STATES 



We now leave the orthonormal boundary states |Q 7 ) and impose periodic boundary 
conditions on the Hamiltonian in Eq. (|l|). We return to the states \Q) n as defined in Eq. (|]), 
where Q is a general m x m matrix, to make an ansatz for the low lying excited states. For 
a translationally invariant system we can define our states to be Bloch states. A reasonable 
ansatz for a Bloch state \Q, k) n defined by a matrix Q and a momentum k is given by 

\Q, k) n = J2T, ^ tr ( A l s n\ ■ ■ ■ A[s j+1 ]QA[ Sj } . . . A[ Sl ] )\s n ... Sl ). (37) 

{%} 3 

This wavefunction can be viewed as the ground state |1) with a disturbance Q introduced at 
some site, and then letting the disturbance run over all sites to form a state with a definite 
momentum. In this way we get a single "particle" excitation. 

As was done for the boundary states in Sec. ^ we now derive expressions for expectation 
values of operators in the states \Q,k) n . The calculations are more tedious and we have 
therefore put the details in App. |^j - [D]. The resulting expressions are in principle similar to 
the ones we obtained for the boundary states, e.g. Eq. fl2Tp. For the norm we find 

(<?', k\Q, fc) = ntr| (Q' <g> 11) £ e ijk l n ' j (TL <g> Q)V j , (38) 

with a similar, but more complicated, expressions for the Hamiltonian and for the z- 
component of total spin. The results can be found in Eq. ( |A1(J| ) and ( |A22| ). We see that 
the general structure of all these matrix elements are that they consist of traces with a 
convolution sum over matrix products inside each trace. For finite length chains, the sums 
in these expressions as well as those in Eq. (|2l[) and (|22D can be expediently calculated by a 
recursive scheme for the case when n is a power of two. These recursive formulas are derived 
in App. [B[ 

One can also calculate the norm and Hamiltonian matrices, G(k, n) and H(k, n), defined 
through the formula 

(Q',k\H op \Q,k) n = nQ'H(k,n)Q (39) 
(Q',k\Q,k) n = nQ'G(k,n)Q, (40) 
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similar to H(ri) and G(n) in Eq. (|T^) and ( p3|) for the boundary states. This time they will 
however depend on k as well as on n. A matrix S^(k,n) representing the z-component of 
total spin can be defined analogously. The principles for calculating these matrices are the 
same as for the boundary states, i.e. one uses the tilde transformation of Eq. (|1^). Due to the 
number of terms in the expression for the expectation values it is numerically cumbersome 
for finite length chains. 

There is however an elegant way to extract the leading behavior of H(k,n) and G(k,n) 
as n — > oo. The details of these calculations can be found in App. |C| and [D[ In this 
section we will only give a brief summary of the method and the results. Let us first 
define the ^-transform (sometimes called a discrete Laplace transform) of a series {a„}^L 
by F(\) = J2^Lo a n£~ nX ■ Let us now denote the sum inside the trace in Eq. (^) by S n , so 
that 

n-l 

S n = e^i n ~ J (l <8>Q)R 

j=0 

We now define a series {S n }^LQ, and take the z-transform of this sequence. By examining 
the analytical structure of the transformed series we are able to extract the leading behavior 
of the sum S n , as n — > oo. In this way we get the asymptotic form of the norm in the limit 
of large n. This procedure is then applied to all sums in the expressions for the matrix 
elements. In App. |C], the ^-transform of a general sum is taken and its large n behavior is 
extracted. In App. [D] we apply the results of App. |C| to the expressions for the expectation 
values derived in App. |A]. This whole procedure finally results in the asymptotic forms 

H(k, n) = n 2 H 2 (k) + nH x {k) + H (k) + 0(z) n (41) 
G(k, n) = nGi(k) + G {k) + 0(z) n , (42) 

with H(k, n) and G(k, n) as defined in Eq. (^) and (|40D. Here z represents the next leading 
eigenvalue of 1 and we find numerically that \z\ ~ 0.8. There are thus very small corrections 
to the asymptotic form. We also find that H 2 and Gi are nonvanishing only when the 
momentum k is zero. The eigenvalue equation which must be solved is 
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H(k,n)Q(k,n) = E(k,n)G(k,n)Q(k,n), 

where Q(k,n) is an m 2 -dimensional vector. For k ^ we thus have 

(nHi(k) + H {k)) Q(k, n) = (nE + A k {n))G {k)Q{k, n), (43) 

where E is the ground state energy per site and Akin) is the excitation energy. E denotes 
the ground state energy per site in the limit n — > oo, and is therefore independent of n. 
Since we are interested in the solutions to Eq. (|43"1) when n — > oo we assume Q and A^ to 
be independent of n and we thus need to solve the simultaneous equations 

H x {k)Q{k) = E G (k)Q(k) (44) 
H (k)Q(k) = A k G (k)Q(k). (45) 

Solving Eq. Q4"5p yields a set of eigenstates Q(k) and eigenvalues A fc for each k. These 
eigenstates have to be simultaneous eigenstates to Eq. flliD with the fc-independent eigenvalue 
Eq. This is in general impossible, unless H\ oc Go, as indeed happens. We thus recover E 
by the proportionality constant 

H x {k) = E G (k). 

The excitation spectrum is then given by the single eigenvalue equation 

H (k)Q(k) = A k G (k)Q(k). (46) 

Similar formulas can be obtained for k = 0. Note that Eq. ( |46D is an eigenvalue equation for 
the excitation spectrum which makes no explicit reference to a ground state. The ground 
state enters however implicitly in the parameters in A[s] on which H(k,n) and G(k,n) 
depends. 

An asymptotic form for the z-component of total spin, S^(k,n), similar to the form for 
H(k,n), containing terms up to order n 2 is also derived in App. 0. Numerically we find 
however, that the only nonvanishing term in S^(k,n), for any momentum k, is the constant 
term S^k). 
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How do we find the orthonormal set of states Q(k) for a particular k from the eigenstate 
equation in Eq. (pp? We can in principle take over the discussion of the boundary states \Q) 
from Sec. |TV|. The only slight problem that enters here is that Go(k) is singular for k ^ 0, that 
is, the nullspace of Go(k) is nonvanishing. In order to find the inverse of Go(k) the nullspace 
must be excluded from the Hilbert space. We do this numerically using singular value 
decomposition. Once this has been done we can simply take over equations (|29|)- (|36|) . In this 
case we identify Ho(k) with H, Go{k) with G and A k with E. We diagonalize Go{k) with V k 
so that V k ;Go(k)Vk — Da{k) is diagonal and define u k = Vk{DG{k))~ 1 ^ 2 . Then we diagonalize 
hk + where h k = u k H (k)u k and s k = u k SQ(k)u k so that w k (h k + es z k )wk = + em k is 
diagonal. We then find 

HoWQ^k) = Ak^G (k)Q^(k), 

where 

{Ql{k)) a ,P = {u k Wk){ a) p) n 

are m x m matrices labeled by 7, orthogonal with respect to Go(k) and simultaneous eigen- 
states of H (k) and S (k). There are less than m 2 eigenvectors Q-y(k) for k 7^ due to the 
nonvanishing nullspace of Go(k). However, probably only a few of the lowest lying energy 
eigenstates Q 1 {k) are reasonable estimates of true excited states. Finally, we can write for 
the orthonormalized states (7, k), defined by the matrices Q^{k) 

| 7 ,/c> = |Q 7 (A;),A;). (47) 
Because states with different values of k are guaranteed to be orthogonal, we find 

(7', k'\j,k) = 6y j7 5 k , jk . 

These represent our "single magnon" states. In the next section we numerically determine 
these states along with their energy and spin expectation values. 
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VI. RESULTS 



We have tested the calculations on the spin-1 Heisenberg chain defined in Eq. (p. All 
computations are done with m = 12, i.e. keeping the twelve states as discussed in Sec. |T|. 
The resultant eight parameter family of trial ground states (Eq. (]?])) was explored. The 
projection matrices A[s] defining the ground state were computed by minimizing the energy 
of the trial ground state. The projection matrices obtained by this variational technique was 
found to agree up to numerical accuracy with the projection operator obtained from similar 
DMRG calculations. The result for the lowest energy state for some (3 is found in Table |. 
The best result known to us for (3 = comes from DMRG calculations in Ref. |||. The exact 
result at the AKLT point (3 = — arctan(l/3) can be found in Ref. || The (3=1 system 
was exactly solved using Bethe ansatz in Ref. |15| The parity operator of Eq. ( |15|) has been 
computed in all cases and it is found that the ground state has parity (— l) n , where n is the 
number of sites. For the string order parameter of Eq. (|T3"D we find g(oo) = —0.3759, whereas 
best estimates are! g(oo) = —0.374325096(2). We find the next leading eigenvalue of 1 to 
be —0.777, giving an asymptotic spin-spin correlation length from Eq. (|14"D of I = 3.963, 
compared to best estimates^ of I = 6.03(1). We believe that the severe truncation of our 
basis to only twelve states has resulted in the asymptotic correlations being quite poor, 
although we have verified that intermediate length spin-spin correlations are consistent with 
more precise calculations^. 

An important issue is whether or not Eq. ([37]) is a good ansatz form for the excitations. 
We have computed the asymptotic forms when n — > oo for the Hamiltonian and norm 
matrices defined in Eq. (|39|) and (flOD as well as for the total spin matrix for different (3 
and momenta k. The orthonormal eigenstates of Eq. (|47| ) are also determined, giving the 
single magnon excitations of our model. The energy and ^-component of total spin for each 
eigenstate are also determined. A particularly interesting point is (3 = 0, the pure Heisenberg 
model, which has been subject to much numerical effort. We find the single particle spectrum 
shown in Fig. [3[ The low-lying triplet branch defines the gap A n = 0.4094, which is very 
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good compared to the most accurately known resultffiO of 0.410502(1). Furthermore, we 
compute the spin wave velocity v = 2.452 to be compared to the calculations in Ref. [7], 
where v = 2.49(1) was obtained. Clearly we reproduce the single-particle triplet excitations 
with high accuracy considering the few number of states in our basis. Our calculation also 
yields a detailed spectrum of lowest lying "single magnon" excitations shown by dotted lines 
in Fig. [3]. Our second lowest energy excitation at k = it is a singlet shown by a dotted line 
in Fig. H with A n (singlet) = 2.348. As a function of k, the second lowest single-particle 
excitation is either a singlet or a spin-2 object, as has also been observed in exact finite 
size calculations^. Parity of each of the elementary excitations is verified by checking the 
relation Eq. (|i~5|) with Q as well as with the matrices A. The boundary to two particle 
excitations at a given value of k is shown in Fig. 0, computed explicitly by minimizing the 
sum of energies of excitations whose pseudomomentum sums to k, and similarly for the three 
particle excitations. These results are shown by the light and dark shaded regions in Fig. |3|. 
The picture fits well with previously obtained results. 

We have similarly computed spectra for various values of /fH^l The result for the gap 
to the lowest lying triplet at k = tt is shown in Fig. [|. Near (3 = 0.6, the excitation spectrum 
at k = 7r crosses zero and becomes negative. Our interpretation of this is that our ground 
state ansatz is deficient, and this shows up as a condensation of elementary excitations. It 
is to be noted that Oitmaa et al.0 also found that numerically the gap appeared to vanish 
rapidly near this value of (3, although they too were unwilling to conclude that this persisted 
in the thermodynamic limit. 

Our calculations are consistent with two possible scenarios of what happens near (3 = 0.6. 
A special value of (3 could exist where the gap closes and signals a new phase. Or, the gap is 
in fact small and persists all the way to (3 = 1 but we do not see it due to our restricted ansatz 
for the ground state. Recent DMRG calculations^ have shown to have similar difficulties 
to estimate the vanishing gap for (3 close to 1. A significant issue appears to be that the 
DMRG fixed point seems to invariably lead to a matrix product ground state that, although 
it succeeds in reproducing ground state energies to high accuracy, cannot strictly give a 
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power law decay of spin correlations. Thus, we find the ground state energy very accurately 
at the Bethe ansatz point (3 = 1 without finding the expected powerlaw decay of correlations. 
The correlation length spectrum is given by the eigenvalues^ of the matrix 1, and it is hard 
to see how this can ever give algebraic correlations. However, intermediate correlations for 
intermediate lengths appear to be well represented in all cases. 

The appendices, App. 0-0, contain the detailed derivations of the results presented in 
Sec. 

VII. CONCLUSIONS 

The present work suggests that the rapid convergence of the DMRG is explained by the 
fact that the states selected are optimally chosen eigenstates of total block spin. Properly 
chosen, these states are highly efficient for building wave functions with a small basis that 
have low total spin for all subblocks. 

Our analysis also proposes that DMRG inherently predicts exponential decay of correla- 
tions. Nevertheless, fully performed DMRG calculations on systems with power law decay 
of correlations seems to agree well with theory. How this is consistent with our calculations 
is currently under study. 

A related topic is the difficulty to describe the vanishing of the gap close to a gapless 
point. However, also "full" DMRG calculations seem to suffer from this problem0. 

APPENDIX A: EXPECTATION VALUES IN THE BLOCH STATES 

In this appendix we will derive expressions for expectation values in the trial Bloch states 
\Q,k) of Eq. © 

n 

\Q, k) = J2 e%km tr ( A \ s n] ■ ■ ■ A[s m+1 ]QA[s m } . . . A[ 8l ] )\s n ... Sl ). 

m=l 

Note that the summation over spins as well as the subscript n, the number of lattice sites, 
are not explicitly written out. 
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1. Calculation of the normalization matrix 

We will derive expressions for expectation values of three types of operators. First we 
calculate the norm (Q f , k\Q, k). Then we show how to obtain the expectation value of total 
spin, St = J2i Si, where Sj = (Sf, Sf, Sf), i.e. the expectation value of the sum of a single 
site operator. Finally we calculate the expectation value of the sum of a two site operator 
like the energy H — J2i$i • The calculations of these three types of expectation values 
differ only in details and not in any fundamental way. For completeness all three cases are 
nevertheless covered in this appendix. 

We begin by calculating the norm of the states \Q,k). Due to the periodic boundary 
condition, states with different k are orthogonal. Using the definition of \Q, k) we have for 
the same value of k, 

(Q', k\Q, k) = E E e- ikm 'e* km tr ( A*[s n ] . . . A*[s m , +1 ]Q'A[s m ,} . . . A*[ Sl ] ) 

m m' {s} 

x tr ( A[s n ] . . . A[s m+1 ]QA[s m ] . . . A[ Sl ] ). (Al) 

We now use the periodic boundary conditions, put w! — 1 and change the summation index 
m - 1 -> m. Using the identities tr ( A ) tr ( B ) = tr ( A ® B ) and (ABC) <g> (DEF) = 
(A <8> D)(B ® E)(C ® F), where the tensor product is defined in the text after Eq. (H), we 
get 

n-l 

(Q', k\Q, k) = n E e lkm tr((A*[s n ] ® A[s n ]){A*[s n ^} A[s n ^}) . . . (A*[s m+2 ] A[s m+2 \) 

m=0 

x(l ® Q)(A*[s m+1 ] ® . . . (Q'* ® ® A[sJ)). 

By defining 

.Rq = 11 ® Q 
L Q = Q* ® 11 

and usine the definition 1 = £ s A*[s] ® A[s] from Eq. (0) we can rewrite this as 
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n-1 

(Q',k\Q,k) =nJ2 e tkm ti ( l n - m - l R Q l m L Q ,i ) 

m=0 

n 

= n £ e tkm ti(L Ql l n - m R Q l m )- ne lkn tr ( L Q , R Q 1 n ) , 

m=0 

where we in the last step have added and subtracted the term m = n and used the cyclicity 
of the trace. Since e lhn = 1 we can now write the norm 

n 

(Q',k\Q,k)=n J2 tr(L Q ,r- m R Q (e tk i) m )- tr(L Q ,R Q l n ) 

m=0 

= ™tr( L Q , (j2(i n - m R Q (e lk i) m ) - RqV^j ) . (A2) 

Let us now introduce the symbol H to represent convolution sums like the one that appears 
inside the trace in Eq. ( |A2| ) . Thus, define the two partition sum S n (x, M,y) by 

E n (x,M,y)= f2x m My n - m , (A3) 

m=0 

where x, M and y are, in our case, square matrices. Later on in this section also three 
partition sums will appear, therefore define 

n n 

Z n (x,M,y,N,z)= J2 E x m 'My m ^Nz n - m \ (A4) 

mi=0 m 2 =mi 

Note that the same symbol, H, is used to represent both two and three partition sums; 
the number of arguments of H determine the number of summation variables. Using this 
definition, the norm can now be written as 

(Q',k\Q,k) = n tr( L QI (E n (i,R Q ,(e ik i))-R Q i n ) )• (A5) 

It is easy to show that Rq and Lq commute, so there is no ambiguity in the order we place 
the Q and the Q' in terms with m = m! in Eq. (|Al|). 

2. Calculation of the total spin 

After finding the norm, we are now interested in the total spin St = SiSj. We thus 
need an expression for the expectation value of the single site operator, 
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(Q',k\(ST)op\Q,k) = ^(Q'^KSiUQ^). 

i=l 

The periodic boundary conditions imply that (Q', k\(Si) op \Q, k) is independent of i so let us 
take i — 1. We then have 

n n 

(Q', k\(S T ) op \Q, k) = nJ2 E e~ ikm 'e lkm tr ( A*[s' n ] . . . AVm'+ilQ^Vm'] . . . AVi] ) 

m=l m'=l 

x tr ( A[s n ] . . . A[s m+ i}QA[s m ] . . . A[si] ) (s' n . . . s\ 

l(Si) op I • • • 

To rewrite this expression using 5 defined in Eq. ( |A3| ) and (|A4|) we split the sums over m 
and ml in three partial sums 

l<m<m' 

Zb= E 

l<m'<m 

Sf= E 

l<m=m' 

Observing the double counting that appear above we see that 

n n 

E E = + S B — £_p. 

m=l m'=l 

Define S^, £b and £p to be the parts of (Q', k\(ST) op \Q, k) with values of m and m! corre- 
sponding to the sums Ha, £b and Hp respectively In a similar way as for the norm we now 
get for the sum A 

n n 

Sa = uJ2 E e^ m - m 'hr(A*is n }...A*[s m/+1 }Q'*...A*ls 2 }A*[s' 1 }) 

m=l rn'=m 

x tr ( A[s n ] . . . A[s m+1 }Q . . . A[s 2 }A[ Sl ] )(s' 1 |(S 1 ) op |s 1 ) 

n n 

= ™E E e lk ^- m 'hi((A*[s n ]0 A[s n ])...(A*[s m , +1 ]® A[s m , +1 ])(Q' ® 1L) 

m=l m'=m 

x(A*[s m ,) ® A[s m ,]) . . . (A*[s m+1 ) ® ® Q) 

x(A*[s m ] ® A[s m ]) . . . GA> 2 ] ® A[s 2 ])(^.*[s'i] ® (^(Si)^*) 

n n 

= n E E e ifc(m " m,) tr ( l n - m ' L Q ,l m '~ m R Q l m - l S ) 

m=l m'=m 
n n 

= n E E e ifc(m - m,) tr ( L Q ,l ml - m R Q l m - x Sl n - m ' ), (A6) 

m=l m'=m 
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where we have used the definition of S from Eq. (0). By changing summation index 
m — 1 — > m and m! — 1 — > m! we get 

n— 1 n— 1 

S A = n £ Y, e ik{m - m ' ] ti ( L Ql l m '- m R Q l m Sl n - l - m ' ) 

m=0 m'=m 

= ntr(L Q ,S ri _ 1 (e- lfc J, i2g, !,£,!)). 
In a similar way we get for the sum B 

n n 

Sn = nJ2 E e ifc ("') tr ( A[s n ] . . . A[s m , +1 ]Q' . . . A[s 2 ]A[s\] )* 

m'=l m=m' 

x tr ( A[s n ] . . . A[s m+1 }Q . . . A[s 2 ]A[ Sl ] )(s / 1 |(S 1 ) op |s 1 ) (A7) 
= ntr(L Q E n ^(l,S,i,R Q ,e ik i) ). (A8) 

It is also possible to show that 

Sb(Q',Q) = (S A (Q,Q'))\ 

The sum F contains the terms that are counted twice in A and B and Sp should therefore 
be subtracted from S A + Sb- We get 



S F = nY tr(l n - m L Q ,R Q l m - l S 

m=l 
n-1 

nY tr(L Q ,R Q l m Sl n - l - m ) 

m=0 

tx(L Q/ RQE n ^(iJ,i)). (A9) 



We now collect the results from Eq. ( |A6| ) , ( |AS| ) and ( |A9[ ). The expectation value of the total 
spin is thus 

(Q', k\(S T ) op \Q, k) = S A (Q', Q) + S B (Q', Q) - S F (Q', Q) 
= n tr ( L Q , x 

[E^ie-^lRg, 1,S,1) 
+ E n _ l (i,S,i,R Q ,e ik i) 

-RgE^lS,!)]). (A10) 
We have here not made use of the fact that Sb can be determined from Sa*- 
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3. Calculation of the energy 



The final operator we need is the energy H = J2i where h^i+x = Si ■ Sj+i. We thus 

have to find an expression for the expectation value of a two site operator. The procedure to 
find it is analogous to how we found the total spin. We use the periodic boundary conditions 
to put 2 = 1. Thus 

n n n 

{Q',k\H op \Q,k) = J2Yl E e- ikm 'e ikm ti(A[s n }...A[s m+1 ]QA[s m ]...A[s 1 ]) 

i=l m=X m'=l 

x tr ( A*[s' n ] . . . A*[s' m , + x}Q'A*[s' m/ } . . . A*[s\] )(s' n . . . s'x\hx, 2 \s n ... Sl ). 

Since the terms with m — 1 and/or m! — 1 in the expression above are special in the sense 
that the matrices Q and Q' mix with the operator /i 1)2 , we this time have to split the sum 
in six partial sums 

£4= E 

2<m<m' 

Zb= E 

2<m'<m 

Zf= E 

2<m=m' 

z c = E 

m=Km' 

Zn= E 

m'=X<m 

Ze= E 

m=m'=X 

We note that 

n n 

E E = + £,B — S.F + Sc + £,D + Eg. 
m=l m'=l 

Analogous to what was done for the total spin, we define Ha, Hb etc. to be the parts of 
(Q', k\H op \Q, k) with values of m and m! corresponding to E^, X# etc. The sum Ha for the 
two site operator is very similar to Sa in Eq. (|A6|) for the single particle operator. We find 

n n 

H A = nY, E e^ m ~ m 'hi((A*[s n ] ® A[s n ]){A*[s n „x) ® . . . ® 

m=2 rn' =m 

x(Q' ® l)(A*[s m /] <g> A[s ro /]) . . . (A*[s m+1 ] ® A[s m+1 ])(IL ® Q) 
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x(A*[s m ] <g> A[s m }) . . . (A*[s' 2 ] ® A[s 2 ])(A*[s'i] ® (s' 2 , sM^K si> 

n n 

J2 e ik{m ~ m ' ) ti( K L Q/ i m '~ m R Q i m - 2 SSi n - m ' ), (All) 



n n 

n 

m,=2 m'—m 



where we have used the hat mapping defined in Eq. ( PH| ) for the Hamiltonian matrix SS. 
By changing summation indices m — 2 — > m and m' — 2 — > m', and using the H notation for 
the sum, we get 



n-2 n-2 

-2-m' 



m=0 m'=m 

= ntr ( L Q ,E n . 2 (e~ lk lR Q , 1,SS,1)). 
In a similar way we get for the sum i7 B 

n n 

#B = n E E e ifc(m - m,) tr ( A[a n ] . . . A[*w +1 ]Q' . . . A[s' 2 ]A[s'i] )* 

m'=2 m=m! 

x tr ( A[s„] . . . A[s m+1 ]Q . . . A[s 2 }A[ Sl } )(s' 2 , s\\h lj2 \s 2 , 8l ) (A12) 
= ntr ( L Q ~ n _ 2 (i,S3,1,Rq, e lk l)). (A13) 

It is also possible to show that 

H B (Q',Q) = (H A (Q,Q'))*. (A14) 

The sum F contains the terms that are counted twice in A and B and Hp should therefore 
be subtracted from H A + Hb- In the same way as we found Sp in Eq. ( |A9| ) we now find 

n 



m=2 

n 



H F = nJ2 tx(l n - m L Q ,R Q l m - 2 SS 

m=2 

tr(L Q/ R Q E n _ 2 (iJS : i)). (A15) 



The sums C, D and E contain terms were the matrix Q and/or Q' mixes with the operator 
h 12 . For C we get 

n 

H c = nJ2 e- lk{m '- x) tr(A*[s n ]...A*[s m>+1 ]Q'*...A*[s' 2 ]A*[s , 1 ]) 

m'=2 

x tr ( A[s n ] . . . A[s 2 ]QA[ Sl ] )(s' 2 , s\\h lt2 \s 2 , 8l ) (A16) 

n 

= nJ2 e~ ik{m '- l) tr ( i n ~ m ' L Q , l m '- 2 SR Q S ) 

m'=2 

= ne- lk tr ( L Q ,E n „ 2 (e- tk i, SR Q S, 1) ). (A17) 
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Hn yields 



H D = n £ e lk ^ tr ( A*[s n ] . . . A*[s' 2 }Q' A*[s\} ) 

x tr( A[s n ] . . . A[s m+1 ]Q. . . A[s 2 ]A[si] )(s' 2 , s'i|/i lj2 |s 2 , si)- 

n 

= n J] e ifc(m - 1} tr ( i n - m R Q l m ~- 2 SL Q ,S ) 

m=2 

= ne~ ik tr ( ~ n _ 2 (V* fc J, l) 5 ), 



(A18) 



(A19) 



where we in the last step used that e %kn = 1. One can also show that 



H D (Q>,Q) = (H C (Q,Q')T. 



The "sum" E is just 



H E = n tr(i n - 2 SR Q L Q ,S 
= n tr ( L,QiSl n ~ 2 SRQ 



(A20) 



We now collect the results from Eq. (TATTD , ( TO ), ( gig) , (|AT7|) , (|SISD and (|A"20|) . For the 
whole Hamiltonian we thus have 



(Q', fc|tf|Q, fc) = H A (Q', Q) + H B (Q', Q) - H F (Q', Q) 

+ H C (Q', Q) + H D {Q', Q) + H E (Q', Q) 
= n tr ( L Q > x 

[s n _ 2 (e- 4fc J,i? Q ,i,^,i) 

+ H n _ 2 (J,55,J,i? Q ,e ifc i) 

— Rq H n _ 2 (i, S^s, i) 

+ e- ifc S n „ 2 (e~ ife i,5 J R Q 5,J) 
+ e- lk S(~ n - 2 (e- lk i,R Q ,i))S 
+ Sl n ~ 2 SR Q \). 



(A21) 



(A22) 



We have here not made use of the relations Hb(Q',Q) = (Ha(Q,Q'))* and Hd(Q',Q) 
(H c (Q,Q')y. 
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Eq. (|A2|) , ( |A10|) and (|A22 ) now contain the desired expectation values, expressed in terms 
of convolution sums. These sums can be expediently calculated using recursive relations, as 
we will show in the next section. 



APPENDIX B: CALCULATING THE PARTITION SUMS RECURSIVELY 

Expectation values between the Bloch states \Q, k) can be divided into partial sums with 
the general forms of two-partition and three-partition sums defined in Eq. ( |A3[ ) - The 
number of terms in the two-partition sum with upper limit n is n + 1 while the number of 
terms in the three-partition sum with upper limit n is ("•+ 1 K n + 2 ) _ Both these sums can be 
calculated recursively with a number of operations of the order log(n). For the two-partition 
sum Eq. (|Al|) we find that the sum with upper limit In can be found from the sum with 
upper limit n by 

E 2n (x, S, y) = x n E n (x, S, y) + E n (x, S, y) y n - x n Sy n 
with the starting sum 

Si (a:, S,y) = xS + Sy. 

We thus get sums where n = 2 J , j integer. Each recursion step requires a constant number 
of additions and multiplications which implies a total computational effort of order log(n). 
The three-partition sum, Eq. (|A4|) , can be done in a similar way. Here the 2n — 2 sum is 
reached from the n — 2 sum by 

E 2n ~2(x, S, y, T, z) = x n S n _ 2 (x, S, y, T, z) + S n _ 2 (^, S, y, T, z) z n 

+ S n _x(x, S, y) S n _i(2/, T, z) 

< 

E 2n _i(x, S,y) = x n S n _i(x, S,y) + E n _ 1 (x,S,y) y n 

y~yi 2 T\ Th 

with similar expressions for S 2n _i(?/,T, z), y 2n and z 2n . Here we start with 
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Eo(x,S,y,T,z) = ST 

Si (a;, S,y) = xS + Sy 
^i(y,T,z) = yT + Tz 

< 

and we get sums with upper summation bound n — 2, with n = 2 J and j an integer. Also 
here the computational effort is of order log(n). In this recursion scheme we also get the 
two-partition sum with upper bound n — 1. 

APPENDIX C: THE POLE EXPANSION 

Although calculating the sums recursively is a nice method for finite size chains, we 
would like to calculate the expectation values in the limit n — > oo. As we will show in this 
section, it is actually possible to do this directly by analyzing the sums' asymptotic form. In 
the next section, App. |D|, we apply the results to the actual sums in the expectation values 
of App. [A]. 

1. Three-partition sums 

In App. |A] expectation values were calculated and expressed in terms of sums. These 
sums are of the general form 

Sn ^ ] (T^) Tx 5 nini -|_ n2 +ri3 ; 

ni,n2,ri3>0 

where x, S and T are m 2 x m 2 matrices and 7 = e tk is a phase factor. We would like to know 
the asymptotic form of S n as n —>■ 00. This form can be found if we take the z-transform 
(also known as discrete Laplace transform) of S n and then analyze the pole structure of the 
transformed sum. Define the ^-transform of the sum 5* n by 

00 

F S [X] ^Y,^^ 
n=0 

We then have 
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^ s w = e E (\ixrs(\xrT(\xr5 n , ni+n2+na 

n=0 ni,ri2,n3 

= E {^lx) ni S{\x) n2 T{\x) n3 

oo oo oo 

= ( E (^Ds( E (^)" 2 )t( E (a^)" 3 ). 

ni=0 ri2=0 ri3=0 

Let us define U as the matrix that diagonalizes x. Let us also define a transformation M 
of a general m 2 x m 2 matrix M by 

M D = U~ l MU. 



Thus x is a diagonal matrix with the eigenvalues of x on the diagonal, while the transfor- 
mation M D of a general matrix M need not be diagonal. We then have 

(oo \ I oo \ I oo \ 

E (A 7 x D ) ni [/-^C/ E (A^)" 2 E (A^)" 3 

ni=0 / \n 2 =0 / \n 3 =0 / 







1 


0^ 




( i 






( 1 


\ 






1 — \~fX\ 




1— A^i 




rpD 


1— \xi 






I 











'''J 




V 


■ / 



where Xj are the eigenvalues of x. In our case x D are the diagonalized 1 and Xj are eigenvalues 
of 1. The largest eigenvalue of I is x x = 1 and the other eigenvalues have absolute values 
less than 0.8. The order of the poles of F s [\] will be different for k = and k ^ 0. We will 
therefore have to treat these two cases separately. We first determine the asymptotic form 
in the k = case. 



a. Pole expansion for zero momentum 

The transform will now have as elements 

{F s W „ = pp, WU _^^_J_. (C1 ) 

Note that we have for simplicity not written out the leading U and the trailing U~ l in the 
above formula. Also in the rest of this article, these U and U~ l will be omitted. Since the 
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largest eigenvalue of 1 is x\ — 1 and the next highest is x<i ~ 0.8, we take as an ansatz for 
the behavior of S n for large n 

S n = An 2 + Bn + C + corrections, 

where the corrections are of order ~ 0.8" and thus very small. We now calculate F [A] 
using this asymptotic form of S n . Call it F A [A] to distinguish it from the original form. 

oo 

F A [A] = A n (An 2 + Bn + C) 

n=0 



AY,n 2 \ n + BY,n\ n + CY,X 



+ (A + B) a 



:i-a) ; 



i-xy 



+ c- 



l-A) 



2A B-3A C -B + A 



(C2) 



(1-A)3 (1-A) 2 (1-A) 

We see that -F A [A] in Eq. (|C2|) has poles at A = 1. F s [\] also has poles at A = 1 and is 
analytical in a neighborhood. We therefore expand F [A] around A = 1 and identify terms. 
This will also justify the asymptotic form we have suggested above. Noting that X\ = 1 we 
define a function g(X) by 



We then note that 



g(\) = (l-X)(l-Xx V ) 



D\-l 



( 



1-A 
1— Ax2 



\ 



(C3) 



(1-A)^[A]=^(A)^(A)T^(A). 

We use the shorthand notation g = g{l) , g' = g'(l) and g" = g"(l)- Combining Eq. (|C2[ ) 
and ( |C3] ) we arrive at the central result of the pole expansion for the three-partition sum 
when k = 0: 
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2A 

-(B-3A) 
2(C - B + A) 



We note that 



lim^ ((1 - A) 3 i^ 
gSgfg 

li^i (^(1-A) 3 F*[A]) 
g'S D gT D g + gS D g'T D g + gS D gT D g' 
lim A^i Qfr(l - A) 3 F*[A]) 
g"S D gT D g + gS D g"T D g + ^ D <?T V 
+ 2gS D g'T D g' + 2g'S D gT D g' + 2g'S D g'T D g. 



( 



9\\) 



\ 



-(l-a=2) 
(1-Ax 2 )2 



-2a 2 (l-a 2 ) 
(1-Ax 2 ) 3 



(C4) 



6. PoZe expansion for nonzero momentum 

We now treat the case when the crystal momentum k 7^ 0. In this case the first matrix 
x is multiplied by a phasefactor 7 = e* fe 7^ 1 and we have the elements 

(F s [x}Y' j = y(s D Y' l (T D y^- — \ — -- — — (C5) 

1 U; Y (1 - A7:*) (1 - Ay,) (1 - A*,-) 1 ' 

We notice that this time there can be no poles of order three at A = 1. Instead we have a 
pole at A = 7 _1 . The asymptotic form now looks like 

S n = Bn + C + 7 n C" + corrections. (C6) 

The new term C will give rise to a term (A — 7~ 1 )~ 1 and to match this term we have to 
expand around A = 7" 1 . There can only be a simple pole at A = 7 _1 so there will not be 
any terms A' or B' (i.e. terms proportional to / ~f n n 2 or j n n). We have 
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F A [A] = ]T xn ( Bn + c + 7 nc ") 



n=0 



B B-C 
+ 



-a 



( A _l)2 (A _!) (At _!) 

B C - B C" 
+ v, rr + 



By defining a function h(X) 



(l-A) (i-a 7 )- 



h(X) = (l-Xx )- 1 

( 1 



(1-Ax 2 ) 



V 



/ 



and using the definition of g(X) we write F S [X] in the following two ways 

1 



F b [\] 



[1 - X) 
1 



-h^l)S u g{X)T u g{X) 



(1 - A 7 ) 

In a similar manner to the k = case we now find 

B = lim A ^((l-A) 2 F s [A] 
= h^gTZg 
-(C-B) = \im x ^U(l-XyF s 



g(X 7 )S v h(X)T»h(X). 



= >yh'(r[)S D gfr D g + h( 1 )S D g'T D g + h( 1 )S D gT D g' 
C = lim^-, ((1-A 7 )F S [A]; 
= gS D h{T l )T D h{T 1 )- 



(C7) 



(C8) 



2. Two-partition sums 

The pole expansion can of course also be done for the two-partition sums defined in 
Eq. (|A3|). We will not go through the details since the calculation is analogous to the 
three-partition case but for completeness only list the results. 

Let us analyze the sum 
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s n = E M m sx n ~ m , 

m=0 

where j,x and S are defined as before. For the case 7 : 
n — > 00 is 

S n = Bn + C + corrections, 

with 

5 = ^ 

< 

-(C-B) = g'S D g + gS D g> 
For the case 7 ^ 1 we instead get the asymptotic form 

S n = C + 7™C" + corrections, 

with 

C = h^)S D g 
C = gS^iT 1 )- 

APPENDIX D: EXPECTATION VALUES USING THE POLE EXPANSION 

In App. ^ we derived expressions for the expectation values of various operators in the 
Bloch states \Q, k). We found that all expectation values were expressed in terms of sums of 
matrix products. In App. |T] we showed that the asymptotic limit of a general sum could be 
calculated. By doing a discrete Laplace transform of the sum and analyzing the analytical 
structure of the transformed sum, we arrived at a closed expression for the asymptotic 
behavior as a sum over just a few matrices. 

In this section we will combine the results of App. [A] and|C| and show how the particular 
sums in the expectation values of App. [A| can be analyzed with the technique of App. 0. By 
doing this we will get rid of the unpleasant sums of App. [A] and replace them with simpler 
expressions describing the asymptotic form of these expectation values in the limit where 
the number of sites goes to infinity. 
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(C9) 

e lh = 1 the asymptotic form as 

(CIO) 

(Gil) 
(C12) 
(C13) 



1. The normalization 



We will begin with the simplest case, the norm as determined in Eq. 

(Q', k\Q } k) = ntr ( Lq,{S° - R Q l n ) ), (Dl) 

with 

In App. y we arrived at two different expressions for the asymptotic forms depending on if 
the momentum k was zero or not. Let us start with k = 0. According to Eq. (|C10| ), the 
sum S„ then has the asymptotic form 



= B G n + C G . 

From Eq. ( |C11|) we directly get 



B G = gR%g 

-(C G -B G ) = g>R°g + gR°g'. 

The last term of Eq. ( pip is no sum and just gives an additional matrix RRg to the asymp- 
totic form of S^. Thus 

(Q', k\Q, k) = n tr ( L Q ,(B G n + C G + Rgg) ) , (D2) 

where k = in this case. Before going on to the case k ^ we will rewrite this formula on a 
more "operator-like" form. This can be done by "pulling out" the matrices Q and Q' from 
the trace. We note that (Q', k\Q,k) in Eq. ( |D2| ) has the form 

(Q', k\Q, k) = tr ( L Q ,M a R Q N a ), 

a 

with M a and N a square matrices on outer product form. By doing a generalization of the 
tilde transformation of Eq. (ITH) we can rewrite this as 
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J2 tr ( L Q ,M a R Q N a ) = J2 Q'*M a N a Q 

a a 

= Q'*feikQV ( J)Q (D3) 

and we find that J2 Q M a N a gives a closed expression for the norm operator, independent of 
Q and Q' (but of course fc-dependent). This transformation can be accomplished by writing 

tr ( L Q ,M a R Q N a ) = (Q 1 * ® mJ-^™^ '(1 ® q^MM^KpMm) 

= q'WqC*) (M/-) ( "' fc) ' (i,m) (jv^)^),^ 
= Q'( i,k ^* Q( n ' p ^ ((M T2341 N T234i ) T3241 )(*' fc )>( n 'P) 
= Q'*A4Jt Q Q } 

where the generalized transpose M Ti ^' h ' 1 is defined in Eq. QT7p. We can thus define a Q' and 
Q independent matrix G(k, n) for k = by 

(Q',k\Q,k) = Q'G(k,n)Q, 



where we determine G(k, n) from Eq. ( |U)2| ) and the generalized tilde transformation Eq. (p3| ). 

Likewise we can derive the expression for G(k,n) for fc ^ 0. This is done in the same 
way by using the formulas Eq. ( |C13|) and (|C12| ). The sum this time has the asymptotic 



form 



S„ — Cq + C'g-, 



where we have assumed n such that e lkn = 1. We now find from Eq. ( p!3 ) that 

C G = gREh(e ik ) 

(D4) 

C'g = h{e-* k )R°g. 

The last term of Eq. ( |D1| ) is independent of k and therefore gives the same contribution, 
Rqg, as before. From Eq. ( p4[ ) and the generalized tilde transformation we can calculate 
G(k, n) also for k ^ 0. 
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2. The Hamiltonian 



Now we calculate the pole expansion of the Hamiltonian in Eq. (|A22j ) . Let us start with 
k ^ this time. We will demonstrate the procedure for the term Ha{Q', Q), just to illustrate 
the three-partition case. For the rest of the terms we will, for completeness, just list the 
results. 

For H A we have from Eq. ( |A22|) 



H A (Q',Q)=ntr(L Q ,Sl 



2 h 



where 



S n -2 — - n _ 2 (e 1 1, Rq, 1, S S ,1 



with the asymptotic form from Eq. ( |C6| ) 



2r< / 



St2 = B A (n-2)+C A + 1 - 2 C 

where we have assumed n such that e tkn = 1. From Eq. ( |C8|) we get 

B A = h{e- lk )R°gS D S D g 
iC A -B A ) = e- lk h\e- lk )R%gS D S D g + h{e- lk )R%g'S D S D g 
+ h{e- ik )R°gS D S D g' 
C A ' = gR°h{e lk )S D S D h{e ik ). 
We now have 

-2/1/ 



H A (Q', Q) =ntr( Lq> (B A (n -2) + C A + A 



We transform this as we did with the norm using Eq. (|D3p to get the i^-operator 

H A (Q',Q) = Q'H A (k,n)Q. 

Note the convention used here. We write the m? x m 2 matrix operator H A , which is in- 
dependent of Q and Q' (but depends on k and n), as H A (k,n) and the (Q', Q)-dependent 
expectation value, which of course also depends on k and n, as H A (Q', Q) . 
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The matrix operator Hb we get from Ha by using Eq. (|A14 ) 



H B (Q',Q) = (H A (Q,Q f )Y 



^Q^M~N a (i ' Mk ' ir Q'^ 

a 

q'* (y^mJjJ) q, 



(D5) 



so that Hs(k,n) = H\(k,n). 

We now do the same procedure for the rest of the sums in Eq. (|A22|) and then later 
also for the case k = 0. The result for the coefficients A, B, C and C in the asymptotic 
expansion of different cases are listed below. Were a coefficient is not present, it is zero. 
Apart from the expression Eq. ( |D5| ) we can also get the asymptotic form of the sum in Hb, 
S^_ 2 = S„_ 2 (J, SS, 1, Rq, e tk l) directly from the pole expansion as 

B B = gR°gS D S D h{e ik ) 
-(Cb-B b ) = e ik gR%gS D S D h'{e ik ) + gR%g'S D S D h{e ik ) 
+ g'R%gS D S D h{e ik ) 
C' B = h{e- ik )R%h{e- ik )S D S D g. 

The sum S^_ 2 = H n _ 2 (J, S D S D , 1) of Hp gives 

B F = gS D S D g 
-{C F -B F ) = g>S D S D g + gS D S D g>. 

The sum S%_ 2 = H n _ 2 (e^ fc i, SR Q S, 1) of H G yields 

C c = h{e- ik )S D R d Q S D g 
C'c = gS D R d Q S D h{e ik ). 

The sum in H Dl S®_ 2 = E n _ 2 (e~ tk l, Rq, 1) can be derived from the relation H D (Q',Q) = 
(Hc{Q, Q'))* in the same way as we did for Hb'- 



H D (k,n) =Hl(k,n), 



or directly from the pole expansion as 
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C D = h(e-* k )R°g 
C' D = gR%h{e* k ). 

Finally, for He of Eq. ( |A2U| ) , which does not contain a sum, we just replace the term i n_2 
by its asymptotic form, g, and then perform the generalized tilde transform of Eq. (p3[). 



The same procedure can be worked out for the Hamiltonian also when k = 0. However 
this time we will have to use the formulas in Eq. flC4[) and (|C11| ) . The technique is analogous 
to the k ^ case and I will not list the results here. 

3. The energy 

Collecting everything together we get for the whole Hamiltonian 

H(k, n) = H A {k, n) + H B {k, n) - H F (k, n) + H c {k, n) + H D (k, n) + H E (k, n) 
and for the energy 

= (Q',k\H op \Q,k) = Q'H(k,n)Q 
Q '' Q[ ' (Q',k\Q,k) Q'G(k,n)Q 

where H(k,n) and G{k,n) are square matrices. This is the result we advertised in Eq. fl39|) , 
(H) and in Eq. © and (0). 



Similar expressions for other expectation values like the total spin in Eq. (|A10 ) can of 
course also be obtained. 
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FIGURES 





n-1 






L tensor product 


1 (3 > t (X) s n 
^ n _l W n 




1 ■ ■ 

L projection 


1 a > 
n 



FIG. 1. A rcnormalization step. First the old block representing n — 1 sites is joined to a single 
site. Tensor products then "glues" the different parts together. Finally there is a projection to a 
new block representing n sites. 
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y j j- p*» r j, 




FIG. 2. Construction of the block states is shown when twelve states are kept in the basis. 
Old representations are on the left and new representations on the right. Each line represents a 
nonzero projection P 1 '' 1 of basis representations. The z-component of total spin is not explicitly 
written out. 
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FIG. 3. The spectrum for f3 = is shown. The lowest single particle triplet is shown as a solid 
line, with the lightly shaded region representing two-particle excitations and the dark region three 
particle excitations. Solid lines define the boundaries to the two and three particle continuum. 
Dotted lines indicate the spectrum of higher energy single-magnon excitations. The spin of these 
dotted excitations are, in order of increasing energy at k = tt: 0,1,2,2,3,1,1,0. 
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A 




FIG. 4. The gap to the lowest triplet excitation at momentum k = tt as a function of is 
shown. 
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TABLES 



(3 Eq exact best numerical 

-\ -0.66666667 § 

-1.40138 - -1.401484038971(4) 

0.6 -2.9184 

1.0 -3.98455 -4 

TABLE I. Ground state energy per site 
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